Case Study: US Craft Beer and Brewery

Introduction

In this case study, we analyzed two data sets provided by our client Budweiser, the Beers data set contains a list of 2410 US craft beers and Breweries data set contains 558 US breweries. The purpose of the code is to find statistical characters of IBU and ABV data, also the correlation between them to address questions from Budweiser.

Repository Structure

Case study Repository in GitHub: https://github.com/ysalame98/Project1 Files: Reers.csv: Source beer data. Breweries.csv: Source breweries data. AlcoholAbuse.csv: US Dept OF HHS https://www.samhsa.gov/data/report/2016-2017-nsduh-state-specific-tables EDAforProject1_R5.RMD: R Markdown file with all R code. EDAforProject1_R5.HTML: HTML file created with Knit.

Analysis:

Data Cleansing and Tidy:

• Identify missing and develop a strategy to resolve
• Identify any duplicates entries
• Since Breweries file has city names look for slight name variations and resolve

Initial Data Profiling

Q1: How many breweries are present in each state?

The number of breweries in 51 states has been calculated and provided in the below table along with a column plot and Tree Plot

# Count breweries in each state
table(BreweriesDF$State)

## 
## AK AL AR AZ CA CO CT DC DE FL GA HI IA ID IL IN KS KY LA MA MD ME MI MN MO 
##  7  3  2 11 39 47  8  1  2 15  7  4  5  5 18 22  3  4  5 23  7  9 32 12  9 
## MS MT NC ND NE NH NJ NM NV NY OH OK OR PA RI SC SD TN TX UT VA VT WA WI WV 
##  2  9 19  1  5  3  3  4  2 16 15  6 29 25  5  4  1  3 28  4 16 10 23 20  1 
## WY 
##  4

2. Merge beer data with the breweries data. Print the first 6 observations and the last six observations to check the merged file.

We merged the two data set into one and named it as ‘mergeDF’, the primary key being used is ‘Brewery_id’ from Beer data set, and ‘Brew_ID’ from Breweries data set. We also changed the two columns’ name for clear understanding. The first and last 6 observations were showed there with head/tail command.

mergeDF <- merge(BeersDF, BreweriesDF, by.x = "Brewery_id", by.y = "Brew_ID")
colnames(mergeDF)[2] = "Beer_name"
colnames(mergeDF)[8] = "Brewery_name"
print("The first 6 observations of the merged file:")

## [1] "The first 6 observations of the merged file:"

head(mergeDF, n=6)

##   Brewery_id     Beer_name Beer_ID   ABV IBU
## 1          1  Get Together    2692 0.045  50
## 2          1 Maggie's Leap    2691 0.049  26
## 3          1    Wall's End    2690 0.048  19
## 4          1       Pumpion    2689 0.060  38
## 5          1    Stronghold    2688 0.060  25
## 6          1   Parapet ESB    2687 0.056  47
##                                 Style Ounces       Brewery_name
## 1                        American IPA     16 NorthGate Brewing 
## 2                  Milk / Sweet Stout     16 NorthGate Brewing 
## 3                   English Brown Ale     16 NorthGate Brewing 
## 4                         Pumpkin Ale     16 NorthGate Brewing 
## 5                     American Porter     16 NorthGate Brewing 
## 6 Extra Special / Strong Bitter (ESB)     16 NorthGate Brewing 
##          City State
## 1 Minneapolis    MN
## 2 Minneapolis    MN
## 3 Minneapolis    MN
## 4 Minneapolis    MN
## 5 Minneapolis    MN
## 6 Minneapolis    MN

print("The last 6 observations of the merged file:")

## [1] "The last 6 observations of the merged file:"

tail(mergeDF, n=6)

##      Brewery_id                 Beer_name Beer_ID   ABV IBU
## 2405        556             Pilsner Ukiah      98 0.055  NA
## 2406        557  Heinnieweisse Weissebier      52 0.049  NA
## 2407        557           Snapperhead IPA      51 0.068  NA
## 2408        557         Moo Thunder Stout      50 0.049  NA
## 2409        557         Porkslap Pale Ale      49 0.043  NA
## 2410        558 Urban Wilderness Pale Ale      30 0.049  NA
##                        Style Ounces                  Brewery_name
## 2405         German Pilsener     12         Ukiah Brewing Company
## 2406              Hefeweizen     12       Butternuts Beer and Ale
## 2407            American IPA     12       Butternuts Beer and Ale
## 2408      Milk / Sweet Stout     12       Butternuts Beer and Ale
## 2409 American Pale Ale (APA)     12       Butternuts Beer and Ale
## 2410        English Pale Ale     12 Sleeping Lady Brewing Company
##               City State
## 2405         Ukiah    CA
## 2406 Garrattsville    NY
## 2407 Garrattsville    NY
## 2408 Garrattsville    NY
## 2409 Garrattsville    NY
## 2410     Anchorage    AK

3. Address the missing values in each column.

The missing values can only be found in column IBU (1005 NA’s) and column ABV (62 NA’s), refer to red bars in the first plot. As we don’t want to delete all the missing values, we started with Option 1 to deal with the missing values, in which the state level median values were calculated for both IBU and ABV and used to replace missing values respectively, but we found that the correlation between the IBU and ABV would be heavily impacted due to this replacement.

Option 1: Simply replacing missing values in IBU/ABV with state median results in significant drop on correlation, from 67% to 49% (Refer to #7)

# Check the missing values in the mergeDF
aggr(mergeDF, prop=FALSE, numbers=TRUE)

marginplot(mergeDF[c("ABV","IBU")], pch = c(20), col = c("darkgray", "red", "blue"))

IBUNADF <-  mergeDF %>% filter(is.na(IBU))
ABVNADF <-  mergeDF %>% filter(is.na(ABV))
BOTHNADF <- mergeDF %>% filter(is.na(ABV) & is.na(IBU))
# Check which state give the most NA's
summary(IBUNADF$State)

##  AK  AL  AR  AZ  CA  CO  CT  DC  DE  FL  GA  HI  IA  ID  IL  IN  KS  KY 
##   8   1   4  23  48 119  21   4   1  21   9   9   5  13  52  48   4   7 
##  LA  MA  MD  ME  MI  MN  MO  MS  MT  NC  ND  NE  NH  NJ  NM  NV  NY  OH 
##   9  31  11  20 124   9  13   0  17  29   0  16   6   0   8   3  28  17 
##  OK  OR  PA  RI  SC  SD  TN  TX  UT  VA  VT  WA  WI  WV  WY 
##   8  38  53   7   9   7   1  41  15   5  10  25  45   0   3

summary(mergeDF$State)

##  AK  AL  AR  AZ  CA  CO  CT  DC  DE  FL  GA  HI  IA  ID  IL  IN  KS  KY 
##  25  10   5  47 183 265  27   8   2  58  16  27  30  30  91 139  23  21 
##  LA  MA  MD  ME  MI  MN  MO  MS  MT  NC  ND  NE  NH  NJ  NM  NV  NY  OH 
##  19  82  21  27 162  55  42  11  40  59   3  25   8   8  14  11  74  49 
##  OK  OR  PA  RI  SC  SD  TN  TX  UT  VA  VT  WA  WI  WV  WY 
##  19 125 100  27  14   7   6 130  26  40  27  68  87   2  15

summary(ABVNADF$State)

## AK AL AR AZ CA CO CT DC DE FL GA HI IA ID IL IN KS KY LA MA MD ME MI MN MO 
##  0  0  0  3  1 15  0  0  1  2  0  0  0  0  0  2  0  1  0  0  0  0 11  0  3 
## MS MT NC ND NE NH NJ NM NV NY OH OK OR PA RI SC SD TN TX UT VA VT WA WI WV 
##  0  1  4  0  4  0  0  1  1  1  0  0  0  3  0  0  0  0  6  0  0  0  0  2  0 
## WY 
##  0

# NADF <- data.frame(State=State, ALL=summary(mergeDF$State), IBUNA=summary(IBUNADF$State), ABVNA=summary(ABVNADF$State))
# NADF$IBUNAperc <- NADF$IBUNA/NADF$ALL
# NADF$ABVNAperc <- NADF$ABVNA/NADF$ALL

# Calculate IBU and ABV median value of each state
IBUmedianState <- mergeDF %>% filter(!is.na(IBU)) %>% group_by(State) %>% dplyr::summarise(IBUmedian=median(IBU))
ABVmedianState <- mergeDF %>% filter(!is.na(ABV)) %>% group_by(State) %>% dplyr::summarise(ABVmedian=median(ABV))

# Replace IBU NA's with state median value
mergeDF_no_IBU_NA <- left_join(mergeDF, IBUmedianState, by = "State")
mergeDF_no_IBU_NA$IBU[is.na(mergeDF_no_IBU_NA$IBU)]  <- as.character(mergeDF_no_IBU_NA$IBUmedian[is.na(mergeDF_no_IBU_NA$IBU)])
aggr(mergeDF_no_IBU_NA, prop=FALSE, numbers=TRUE) # Still 7 observations IBU NA, because there is no any IBU value available for state 

# Replace 7 IBU NA's in SD with grand median value 35.00
mergeDF_no_IBU_NA$IBU = str_replace_na(mergeDF_no_IBU_NA$IBU, replacement =  "35")
mergeDF_no_IBU_NA$IBUmedian = str_replace_na(mergeDF_no_IBU_NA$IBUmedian, replacement =  "35")

# Replace ABV NA's with state median value
mergeDF_noNA <- left_join(mergeDF_no_IBU_NA, ABVmedianState, by = "State")
mergeDF_noNA$ABV[is.na(mergeDF_noNA$ABV)] <- as.character(mergeDF_noNA$ABVmedian[is.na(mergeDF_noNA$ABV)])
aggr(mergeDF_noNA, prop=FALSE, numbers=TRUE)

mergeDF_noNA$IBU <- as.numeric(mergeDF_noNA$IBU)
mergeDF_noNA$ABV <- as.numeric(mergeDF_noNA$ABV)

# Plot
IBU.plotly <- mergeDF_noNA %>% ggplot(aes(x=IBU))+geom_bar()+ggtitle("IBU Distribution") +
  theme(axis.text.x = element_text(size = 6, angle = 90))
ggplotly(IBU.plotly)

ABV.plotly <- mergeDF_noNA %>% ggplot(aes(x=ABV))+geom_bar()+ggtitle("ABV Distribution") +
  theme(axis.text.x = element_text(size = 6, angle = 90))
ggplotly(ABV.plotly)

# Statistic ABV, IBU
summary(mergeDF_noNA)

##    Brewery_id                     Beer_name       Beer_ID      
##  Min.   :  1.0   Nonstop Hef Hop       :  12   Min.   :   1.0  
##  1st Qu.: 94.0   Dale's Pale Ale       :   6   1st Qu.: 808.2  
##  Median :206.0   Oktoberfest           :   6   Median :1453.5  
##  Mean   :232.7   Longboard Island Lager:   4   Mean   :1431.1  
##  3rd Qu.:367.0   1327 Pod's ESB        :   3   3rd Qu.:2075.8  
##  Max.   :558.0   Boston Lager          :   3   Max.   :2692.0  
##                  (Other)               :2376                   
##       ABV               IBU                                    Style     
##  Min.   :0.00100   Min.   :  4.00   American IPA                  : 424  
##  1st Qu.:0.05000   1st Qu.: 26.00   American Pale Ale (APA)       : 245  
##  Median :0.05700   Median : 35.00   American Amber / Red Ale      : 133  
##  Mean   :0.05973   Mean   : 39.81   American Blonde Ale           : 108  
##  3rd Qu.:0.06700   3rd Qu.: 47.00   American Double / Imperial IPA: 105  
##  Max.   :0.12800   Max.   :138.00   American Pale Wheat Ale       :  97  
##                                     (Other)                       :1298  
##      Ounces                          Brewery_name      City          
##  Min.   : 8.40   Brewery Vivant            :  62   Length:2410       
##  1st Qu.:12.00   Oskar Blues Brewery       :  46   Class :character  
##  Median :12.00   Sun King Brewing Company  :  38   Mode  :character  
##  Mean   :13.59   Cigar City Brewing Company:  25                     
##  3rd Qu.:16.00   Sixpoint Craft Ales       :  24                     
##  Max.   :32.00   Hopworks Urban Brewery    :  23                     
##                  (Other)                   :2192                     
##      State       IBUmedian           ABVmedian      
##  CO     : 265   Length:2410        Min.   :0.04000  
##  CA     : 183   Class :character   1st Qu.:0.05500  
##  MI     : 162   Mode  :character   Median :0.05700  
##  IN     : 139                      Mean   :0.05677  
##  TX     : 130                      3rd Qu.:0.05800  
##  OR     : 125                      Max.   :0.06250  
##  (Other):1406

Option 2: State level median values were only used to replace ABV missing values (62 NA’s in ABV), then we predicted the IBU values (1005 NA’s) with a Simple Linear Regression model base on the exploratory variable ABV, this way we keep the correlation between IBU/ABV no big change. This is the option that we used to get the data frame ‘tidyDF’, which is used to address the following questions.

# Calculate ABV median value of each state
ABVmedianState <- mergeDF %>% filter(!is.na(ABV)) %>% group_by(State) %>%
  dplyr::summarise(ABVmedian=median(ABV))
# Replace 62 NA's in ABV with state level median
mergeDF_no_ABV_NA <- left_join(mergeDF, ABVmedianState, by = "State")

mergeDF_no_ABV_NA$ABV[is.na(mergeDF_no_ABV_NA$ABV)] <-
  mergeDF_no_ABV_NA$ABVmedian[is.na(mergeDF_no_ABV_NA$ABV)]

aggr(mergeDF_no_ABV_NA, prop=FALSE, numbers=TRUE)

# Build simple regression model with 1405 pairs of IBU/ABV to get the correlation coef.
mergeDF_tidy <- mergeDF %>% filter(!is.na(IBU) & !is.na(ABV))
mod <- lm(mergeDF_tidy$IBU ~ mergeDF_tidy$ABV)
summary(mod)

## 
## Call:
## lm(formula = mergeDF_tidy$IBU ~ mergeDF_tidy$ABV)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -78.849 -11.977  -0.721  13.997  93.458 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       -34.099      2.326  -14.66   <2e-16 ***
## mergeDF_tidy$ABV 1282.037     37.860   33.86   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 19.26 on 1403 degrees of freedom
## Multiple R-squared:  0.4497, Adjusted R-squared:  0.4493 
## F-statistic:  1147 on 1 and 1403 DF,  p-value: < 2.2e-16

Intercept <- mod$coefficients[1]
Slope <- mod$coefficients[2]

# Add new column predicted_IBU base on correlation coef.
# Change negative predicted_IBU to 0, IBU is measured on a scale from 0 to infinite
mergeDF_no_ABV_NA$predict_IBU <- pmax(Intercept + Slope*mergeDF_no_ABV_NA$ABV, 0)
summary(mergeDF_no_ABV_NA$IBU)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##    4.00   21.00   35.00   42.71   64.00  138.00    1005

# Replace 1005 NA's in IBU with predicted_IBU
mergeDF_no_ABV_NA$IBU[is.na(mergeDF_no_ABV_NA$IBU)] <-
  mergeDF_no_ABV_NA$predict_IBU[is.na(mergeDF_no_ABV_NA$IBU)]
summary(mergeDF_no_ABV_NA$IBU)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00   26.00   36.41   42.49   55.64  138.00

tidyDF <- mergeDF_no_ABV_NA
aggr(tidyDF, prop=FALSE, numbers=TRUE)

# makeCodebook(tidyDF)
# makeDataReport(tidyDF)

4. Compute the median alcohol content(ABV) and international bitterness unit (IBU) for each state. Plot a bar chart to compare.

Median ABV and IBU data have been sorted and plotted using bar char with flipped coordinates. Due to large # on bars, this makes the charts hard to read. Line chart with points will be a choice.

# plot IBU median distribution using IBUmedianState data frame

## Calculate IBU and ABV median value of each state
IBUmedianState <- tidyDF %>% filter(!is.na(IBU)) %>% group_by(State) %>% dplyr::summarise(IBUmedian=median(IBU))
ABVmedianState <- tidyDF %>% filter(!is.na(ABV)) %>% group_by(State) %>% dplyr::summarise(ABVmedian=median(ABV))
## Plot
IBUmedianState$State <- factor(IBUmedianState$State, level = IBUmedianState$State[order(IBUmedianState$IBUmedian)])
IBUmedianState %>% 
  ggplot(aes(x=State, y=IBUmedian, fill=State)) + 
  geom_bar(stat = "identity") + 
  ggtitle("IBU Median By States")+
#   geom_text(aes(label=IBUmedian), hjust=-0.5) +
#   theme(axis.text.y = element_text(size = 4)) +
#   coord_flip()
# ggplotly()

geom_text(mapping = aes(label = sprintf("%0.2f", round(IBUmedian,digits = 2)), size=3), check_overlap = TRUE) +
  theme(axis.text = element_text(size =7) ,legend.position = "none") +
  coord_flip()

ggplotly()

# plot ABV median distribution using ABVmedianState data frame
ABVmedianState$State <- factor(ABVmedianState$State, level = ABVmedianState$State[order(ABVmedianState$ABVmedian)])
ABVmedianState %>% ggplot(aes(x=State, y=ABVmedian, fill=State)) + geom_bar(stat = "identity") + ggtitle("ABV Median By States")+
  geom_text(aes(label=ABVmedian), hjust=-0.5) +
  theme(axis.text.y = element_text(size = 6)) +
  coord_flip()

ggplotly()

# geom_text(mapping = aes(label = sprintf("%0.3f", round(ABVmedian,digits = 4)), hjust=-0.5,size=3), check_overlap = TRUE) +
#   theme(axis.text = element_text(size =7) ,legend.position = "none") +
#   coord_flip()
# ggplotly()
# plot scatter with 51 pairs of median values of ABV and IBU
ABV_IBU_mergeDF <- merge(ABVmedianState, IBUmedianState, by = "State" )
ABV_IBU_mergeDF$State <- factor(ABV_IBU_mergeDF$State, droplevels(ABV_IBU_mergeDF$State))
ABV_IBU_mergeDF %>% ggplot(aes(x=IBUmedian, y=ABVmedian)) + geom_point()

5. Which state has the maximum alcoholic (ABV) beer? Which state has the most bitter (IBU) beer?

Maximum ABV and IBU data have been sorted and plotted.

IBUmaxState <- tidyDF %>% filter(!is.na(IBU)) %>% group_by(State) %>% dplyr::summarise(IBUmax=max(IBU))
ABVmaxState <- tidyDF %>% filter(!is.na(ABV)) %>% group_by(State) %>% dplyr::summarise(ABVmax=max(ABV))

# plot IBU max distribution using IBUmedianState data frame
IBUmaxState$State <- factor(IBUmaxState$State, level = IBUmaxState$State[order(IBUmaxState$IBUmax)])
IBUmaxState %>% ggplot(aes(x=State, y=IBUmax, fill=State)) + geom_bar(stat = "identity") + ggtitle("Max IBU By State")+
#   geom_text(aes(label=IBUmax), hjust=-0.5) +
#   theme(axis.text.y = element_text(size = 6) ,legend.position = "none") +
#   coord_flip()
# ggplotly()

geom_text(mapping = aes(label = sprintf("%0.2f", round(IBUmax,digits = 2)), size=3), check_overlap = TRUE) +
  theme(axis.text = element_text(size =7) ,legend.position = "none") +
  ylab("Max IBU By State") +
  coord_flip()

ggplotly()

# plot ABV max distribution using ABVmaxState data frame
ABVmaxState$State <- factor(ABVmaxState$State, level = ABVmaxState$State[order(ABVmaxState$ABVmax)])
ABVmaxState %>% ggplot(aes(x=State, y=ABVmax, fill=State)) + geom_bar(stat = "identity") + ggtitle("ABV max By State")+
#   geom_text(aes(label=ABVmax), hjust=-0.5) +
#   theme(axis.text.y = element_text(size = 6)) +
#   coord_flip()
# ggplotly()

geom_text(mapping = aes(label = sprintf("%0.2f", round(ABVmax,digits = 2)), size=3), check_overlap = TRUE) +
  theme(axis.text = element_text(size =7) ,legend.position = "none") +
  ylab("Max ABV By State") +
  coord_flip()

ggplotly()

6. Comment on the summary statistics and distribution of the ABV variable.

Statistics of ABV variable can be found with summary function, box plot, bar chart and scatter plots. Both ABV and IBU with ABV showing some outliers in the 12oz cans category and IBU exhibits some outliers also in 12 OZ IPA.

summary(tidyDF)

##    Brewery_id                     Beer_name       Beer_ID      
##  Min.   :  1.0   Nonstop Hef Hop       :  12   Min.   :   1.0  
##  1st Qu.: 94.0   Dale's Pale Ale       :   6   1st Qu.: 808.2  
##  Median :206.0   Oktoberfest           :   6   Median :1453.5  
##  Mean   :232.7   Longboard Island Lager:   4   Mean   :1431.1  
##  3rd Qu.:367.0   1327 Pod's ESB        :   3   3rd Qu.:2075.8  
##  Max.   :558.0   Boston Lager          :   3   Max.   :2692.0  
##                  (Other)               :2376                   
##       ABV               IBU                                    Style     
##  Min.   :0.00100   Min.   :  0.00   American IPA                  : 424  
##  1st Qu.:0.05000   1st Qu.: 26.00   American Pale Ale (APA)       : 245  
##  Median :0.05700   Median : 36.41   American Amber / Red Ale      : 133  
##  Mean   :0.05973   Mean   : 42.49   American Blonde Ale           : 108  
##  3rd Qu.:0.06700   3rd Qu.: 55.64   American Double / Imperial IPA: 105  
##  Max.   :0.12800   Max.   :138.00   American Pale Wheat Ale       :  97  
##                                     (Other)                       :1298  
##      Ounces                          Brewery_name      City          
##  Min.   : 8.40   Brewery Vivant            :  62   Length:2410       
##  1st Qu.:12.00   Oskar Blues Brewery       :  46   Class :character  
##  Median :12.00   Sun King Brewing Company  :  38   Mode  :character  
##  Mean   :13.59   Cigar City Brewing Company:  25                     
##  3rd Qu.:16.00   Sixpoint Craft Ales       :  24                     
##  Max.   :32.00   Hopworks Urban Brewery    :  23                     
##                  (Other)                   :2192                     
##      State        ABVmedian        predict_IBU    
##  CO     : 265   Min.   :0.04000   Min.   :  0.00  
##  CA     : 183   1st Qu.:0.05500   1st Qu.: 30.00  
##  MI     : 162   Median :0.05700   Median : 38.98  
##  IN     : 139   Mean   :0.05677   Mean   : 42.49  
##  TX     : 130   3rd Qu.:0.05800   3rd Qu.: 51.80  
##  OR     : 125   Max.   :0.06250   Max.   :130.00  
##  (Other):1406

# Create new column for IPA/Ale breakdown
tidyDF$IPA_Ale <- ifelse(str_detect(tidyDF$Style,regex("IPA", ignore_case = TRUE)),"IPA",ifelse(str_detect(tidyDF$Style,regex("Ale",ignore_case = TRUE)),"Ale","Other"))

# Visualization with ABV data
##   Bloxplot
tidyDF %>% ggplot(aes(y=ABV, x=IPA_Ale, fill=as.factor(Ounces))) + geom_boxplot() +
  xlab("Ale Types (IPA or ALE)") + ylab("ABV") + ggtitle("Boxplot of ABV vs Ale Type")

tidyDF %>% ggplot(aes(y=ABV, x=as.factor(Ounces), fill=IPA_Ale)) + geom_boxplot() +
  xlab("Ounces") + ylab("ABV") + ggtitle("Boxplot of ABV vs Ounces")

##   Histogram
tidyDF %>% ggplot(aes(x=ABV)) + geom_histogram(fill="blue", color="black")

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

##   Barchart
tidyDF %>% ggplot(aes(x=ABV, fill=IPA_Ale)) + geom_bar(color="black") + xlim(c(0.025,0.100)) +
  xlab("Alcohol By Volume (ABV)") + ylab("Beer Name Count") + ggtitle("Ale/IPA Weight vs ABV")

## Warning: Removed 5 rows containing non-finite values (stat_count).

## Warning: Removed 1 rows containing missing values (geom_bar).

tidyDF %>% ggplot(aes(x=ABV, fill=IPA_Ale)) + geom_bar(color="black", position = "fill") + xlim(c(0.025,0.100)) +
  xlab("Alcohol By Volume (ABV)") + ylab("Beer Name Count%") + ggtitle("Ale/IPA Weight% vs ABV")

## Warning: Removed 5 rows containing non-finite values (stat_count).

## Warning: Removed 1 rows containing missing values (geom_bar).

tidyDF %>% ggplot(aes(x=ABV, fill=as.factor(Ounces))) + geom_bar(color="black") + xlim(c(0.025,0.100)) +
  xlab("Alcohol By Volume (ABV)") + ylab("Beer Name Count") + ggtitle("Ounce Weight vs ABV")

## Warning: Removed 5 rows containing non-finite values (stat_count).

## Warning: position_stack requires non-overlapping x intervals

## Warning: Removed 3 rows containing missing values (geom_bar).

# Visualization with IBU data
##   Bloxplot
tidyDF %>% ggplot(aes(y=IBU, x=IPA_Ale, fill=as.factor(Ounces))) + geom_boxplot() +
  xlab("Ale Types (IPA or ALE)") + ylab("IBU") + ggtitle("Boxplot of IBU vs Ale Type")

tidyDF %>% ggplot(aes(y=IBU, x=as.factor(Ounces), fill=IPA_Ale)) + geom_boxplot() +
  xlab("Ounces") + ylab("IBU") + ggtitle("Boxplot of IBU vs Ounces")

##   Barchart
tidyDF %>% ggplot(aes(x=IBU, group = interaction(IBU, IPA_Ale), fill=IPA_Ale)) + geom_bar(color="black") +
  xlab("International Bitterness Units (IBU)") + ylab("Beer Name Count") + ggtitle("Ale/IPA Weight vs IBU")

## Warning: position_stack requires non-overlapping x intervals

tidyDF %>% ggplot(aes(x=IBU, group = interaction(IBU, IPA_Ale), fill=IPA_Ale)) + geom_bar(color="black", position = "fill") +
  xlab("International Bitterness Units (IBU)") + ylab("Beer Name Count%") + ggtitle("Ale/IPA Weight% vs IBU")

## Warning: position_stack requires non-overlapping x intervals

tidyDF %>% ggplot(aes(x=IBU, group = interaction(IBU, IPA_Ale), fill=as.factor(Ounces))) + geom_bar(color="black") +
  xlab("International Bitterness Units (IBU)") + ylab("Beer Name Count") + ggtitle("Ounce Weight vs IBU")

## Warning: position_stack requires non-overlapping x intervals

# Scatter plot ABV vs. IBU
tidyDF %>% ggplot(aes(x=IBU, y=ABV)) + geom_point(shape=0) + ggtitle("ABV vs IBU") + geom_smooth(method = lm)

tidyDF %>% ggplot(aes(x=IBU, y=ABV, color=IPA_Ale)) + geom_point() + ggtitle("ABV vs IBU") + geom_smooth(method = lm)

tidyDF %>% ggplot(aes(x=IBU, y=ABV, color=IPA_Ale)) + geom_point() + ggtitle("ABV vs IBU") + facet_wrap(~State)

tidyDF %>% ggplot(aes(x=IBU, y=ABV, color=IPA_Ale)) + geom_point() + ggtitle("ABV vs IBU") + facet_wrap(~as.factor(Ounces))

7. Is there an apparent relationship between the bitterness of the beer and its alcoholic content? Draw a scatter plot. Make your best judgment of a relationship and EXPLAIN your answer.

Scatter plots with linear model smooth function are provided, the correlation between IBU and ABV is 0.757878, with a p-value < 2e-16, Multiple R-squared is 0.5744, meaning that 57% changes in exploratory variable ABV can explain the changes in response variable IBU.

# Scatter plot and Pearson's correlation with original data
mergeDF %>% filter(!is.na(IBU) & !is.na(ABV)) %>%
  ggplot(aes(x=ABV, y=IBU)) + 
  geom_point(aes(color=State), size=1, position="jitter") + 
  geom_smooth(method = lm) + 
  theme_wsj() +
  xlab("ABV") + ylab("IBU") +
  ggtitle("Relationship between IBU and ABV")

ggplotly()

## Warning: plotly.js does not (yet) support horizontal legend items 
## You can track progress here: 
## https://github.com/plotly/plotly.js/issues/53

mergeDF_tidy <- mergeDF %>% filter(!is.na(IBU) & !is.na(ABV))
cor.test(mergeDF_tidy$IBU, mergeDF_tidy$ABV)

## 
##  Pearson's product-moment correlation
## 
## data:  mergeDF_tidy$IBU and mergeDF_tidy$ABV
## t = 33.863, df = 1403, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.6407982 0.6984238
## sample estimates:
##       cor 
## 0.6706215

mod <- lm(mergeDF_tidy$IBU ~mergeDF_tidy$ABV)
summary(mod)

## 
## Call:
## lm(formula = mergeDF_tidy$IBU ~ mergeDF_tidy$ABV)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -78.849 -11.977  -0.721  13.997  93.458 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       -34.099      2.326  -14.66   <2e-16 ***
## mergeDF_tidy$ABV 1282.037     37.860   33.86   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 19.26 on 1403 degrees of freedom
## Multiple R-squared:  0.4497, Adjusted R-squared:  0.4493 
## F-statistic:  1147 on 1 and 1403 DF,  p-value: < 2.2e-16

# Scatter plot and Pearson's correlation with tidy data (no missing value)
tidyDF %>% ggplot(aes(x=ABV, y=IBU)) + 
  geom_point(aes(color=State), size=1, position="jitter") + 
  geom_smooth(method = lm) + 
  ggtitle("Relationship between IBU and ABV")

ggplotly()

cor.test(tidyDF$IBU, tidyDF$ABV)

## 
##  Pearson's product-moment correlation
## 
## data:  tidyDF$IBU and tidyDF$ABV
## t = 57.005, df = 2408, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.7403535 0.7743731
## sample estimates:
##      cor 
## 0.757878

mod1 <- lm(tidyDF$IBU ~ tidyDF$ABV)
summary(mod1)

## 
## Call:
## lm(formula = tidyDF$IBU ~ tidyDF$ABV)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -78.718  -4.712  -0.035   2.762  93.442 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -33.818      1.372  -24.65   <2e-16 ***
## tidyDF$ABV  1277.567     22.411   57.01   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 14.72 on 2408 degrees of freedom
## Multiple R-squared:  0.5744, Adjusted R-squared:  0.5742 
## F-statistic:  3250 on 1 and 2408 DF,  p-value: < 2.2e-16

#Correlations
Corl <- inspect_cor(mergeDF)
show_plot (Corl) + ggtitle("Correlation Between All Attributes")

# Sctter Plot
mergeDF %>% ggplot(aes(x=IBU, y=ABV)) + geom_point(aes(color=State), size=1, position="jitter") + geom_smooth()

## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'

## Warning: Removed 1005 rows containing non-finite values (stat_smooth).

## Warning: Removed 1005 rows containing missing values (geom_point).

# Plot GGpairs
tidyDF %>%
select(ABV, IBU, IPA_Ale) %>%
ggpairs(mapping = aes(color = IPA_Ale))

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

8. Budweiser would also like to investigate the difference with respect to IBU and ABV between IPAs (India Pale Ales) and other types of Ale (any beer with “Ale” in its name other than IPA). You decide to use KNN clustering to investigate this relationship. Provide statistical evidence one way or the other. You can of course assume your audience is comfortable with percentages … KNN is very easy to understand.

A category variable IPA_Ale is created to split all the 2410 beers into three groups: Ale (1007), IPA (571) and Other (832), we only use the first two groups for analysis in this part. With this KNN classification model, we get an estimated accuracy of 0.7893 with k=5.

table(tidyDF$IPA_Ale) # Ale:963, IPA:571, Other:876

## 
##   Ale   IPA Other 
##  1007   571   832

# Check the Ale: 2nd time with regex(): # Ale:1007, IPA:571, Other:832

tidyDF.knn <- tidyDF %>% filter(IPA_Ale == "IPA" | IPA_Ale == "Ale")

## 
set.seed(1234)
splitPerc = .75
iterations = 50
numks = 30

masterAcc = matrix(nrow = iterations, ncol = numks)
  
for(j in 1:iterations)
{
  accs = data.frame(accuracy = numeric(30), k = numeric(30))
  trainIndices = sample(1:dim(tidyDF.knn)[1],round(splitPerc * dim(tidyDF.knn)[1]))
  train = tidyDF.knn[trainIndices,]
  test = tidyDF.knn[-trainIndices,]
  for(i in 1:numks)
    {
    classifications = knn(train[,c(4,5)],test[,c(4,5)],train$IPA_Ale, prob = TRUE, k = i)
    table(classifications,test$IPA_Ale)
    CM = confusionMatrix(table(classifications,test$IPA_Ale))
    masterAcc[j,i] = CM$overall[1]
  }
}

MeanAcc = colMeans(masterAcc)
maxMeanAcc = max(MeanAcc) # k=?
plot(seq(1,numks,1),MeanAcc, type = "l")

classifications5 = knn(train[,c(4,5)],test[,c(4,5)],train$IPA_Ale, prob = TRUE, k = 5)
table(classifications5,test$IPA_Ale)

##                 
## classifications5 Ale IPA
##              Ale 226  46
##              IPA  31  91

CM5 = confusionMatrix(table(classifications5,test$IPA_Ale))
CM5

## Confusion Matrix and Statistics
## 
##                 
## classifications5 Ale IPA
##              Ale 226  46
##              IPA  31  91
##                                           
##                Accuracy : 0.8046          
##                  95% CI : (0.7619, 0.8426)
##     No Information Rate : 0.6523          
##     P-Value [Acc > NIR] : 2.206e-11       
##                                           
##                   Kappa : 0.5579          
##                                           
##  Mcnemar's Test P-Value : 0.1106          
##                                           
##             Sensitivity : 0.8794          
##             Specificity : 0.6642          
##          Pos Pred Value : 0.8309          
##          Neg Pred Value : 0.7459          
##              Prevalence : 0.6523          
##          Detection Rate : 0.5736          
##    Detection Prevalence : 0.6904          
##       Balanced Accuracy : 0.7718          
##                                           
##        'Positive' Class : Ale             
## 

set.seed(1234)
splitPerc = .80
iterations = 50
numks = 30

masterAcc = matrix(nrow = iterations, ncol = numks)
  
for(j in 1:iterations)
{
  accs = data.frame(accuracy = numeric(30), k = numeric(30))
  trainIndices = sample(1:dim(tidyDF.knn)[1],round(splitPerc * dim(tidyDF.knn)[1]))
  train = tidyDF.knn[trainIndices,]
  test = tidyDF.knn[-trainIndices,]
  for(i in 1:numks)
    {
    classifications = knn(train[,c(4,5)],test[,c(4,5)],train$IPA_Ale, prob = TRUE, k = i)
    table(classifications,test$IPA_Ale)
    CM = confusionMatrix(table(classifications,test$IPA_Ale))
    masterAcc[j,i] = CM$overall[1]
  }
}

MeanAcc = colMeans(masterAcc)
maxMeanAcc = max(MeanAcc) # k=?
plot(seq(1,numks,1),MeanAcc, type = "l")

classifications5 = knn(train[,c(4,5)],test[,c(4,5)],train$IPA_Ale, prob = TRUE, k = 6)
table(classifications5,test$IPA_Ale)

##                 
## classifications5 Ale IPA
##              Ale 182  40
##              IPA  22  72

CM5 = confusionMatrix(table(classifications5,test$IPA_Ale))
CM5

## Confusion Matrix and Statistics
## 
##                 
## classifications5 Ale IPA
##              Ale 182  40
##              IPA  22  72
##                                           
##                Accuracy : 0.8038          
##                  95% CI : (0.7557, 0.8461)
##     No Information Rate : 0.6456          
##     P-Value [Acc > NIR] : 5.52e-10        
##                                           
##                   Kappa : 0.5551          
##                                           
##  Mcnemar's Test P-Value : 0.03085         
##                                           
##             Sensitivity : 0.8922          
##             Specificity : 0.6429          
##          Pos Pred Value : 0.8198          
##          Neg Pred Value : 0.7660          
##              Prevalence : 0.6456          
##          Detection Rate : 0.5759          
##    Detection Prevalence : 0.7025          
##       Balanced Accuracy : 0.7675          
##                                           
##        'Positive' Class : Ale             
## 

9. Knock their socks off! Find one other useful inference from the data that you feel Budweiser may be able to find value in. You must convince them why it is important and back up your conviction with appropriate statistical evidence.

We identified an external data source from the Dept of HHS which shows the level of alcohol abuse by various age group. We hypothesized that there may be some correlation between the abundant availability of alcohol, as represented by number of breweries and the number of cases of alcohol abuse.

Using a map chart we plotted the continuous variables and sure enough there visual evidence that supports our hypothesis. We also conducted a correlation analysis, and we obtained a positive linear relationship

  lookup = data.frame(abb=state.abb, state= state.name )
  #View(lookup)
  # add new column to tidyDF with lower case state codes to user to merge with lookup DF
  abb = unlist(str_split(str_replace_all(tidyDF$State,"[^[:alnum:] ]", ""), boundary("word")))
  tidyDF$abb <- as.factor(abb)
  #View(tidyDF)
  
  # create a new data set with BrewriesDF and lookup data
  BrewriesDF1 = merge(tidyDF,lookup,"abb")
  #View(BrewriesDF1)
  # count # of sate in Acu2 data frame
  BrewriesMPD= count(BrewriesDF1,"state")
  #View(BrewriesMPD)
  # Change the name of the count column to Brewreies Count
  colnames(BrewriesMPD)[2] = "BrewriesCnt"
  #change the names of the states to lower case.  
  BrewriesMPD$region <- tolower(BrewriesMPD$state)
  # get the longi/lat of the states in order to draw the map
  states <- map_data("state")
  #View(states)
  #merge the states long/lat data with AcuMap data
  map.df <- left_join(states,BrewriesMPD, by = "region")
  head(BrewriesMPD[order(-BrewriesMPD$BrewriesCnt),], n=11)

##            state BrewriesCnt        region
## 6       Colorado         265      colorado
## 5     California         183    california
## 22      Michigan         162      michigan
## 14       Indiana         139       indiana
## 43         Texas         130         texas
## 37        Oregon         125        oregon
## 38  Pennsylvania         100  pennsylvania
## 13      Illinois          91      illinois
## 49     Wisconsin          87     wisconsin
## 21 Massachusetts          82 massachusetts
## 32      New York          74      new york

  # Breweries Top10 cities latitude and longitude
  DenverBrew <- tibble(long = c(-104.991531), lat = c(39.742043), names = c("CO\n265"))
  LosAngelesBrew <- tibble(long = c(-118.243683), lat = c(34.052235), names = c("CA\n183"))
  TraverseCityBrew <- tibble(long = c(-85.6206), lat = c(44.7631), names = c("MI\n162"))
  AustinBrew <- tibble(long = c(-97.743057), lat = c(30.267153), names = c("TX\n130"))  
  PhiladelphiaBrew <- tibble(long = c(-75.165222), lat = c(39.952583), names = c("PA\n100"))
  ChicagoBrew <- tibble(long = c(-87.629799), lat = c(41.878113), names = c("IL\n91"))
  BostonBrew <- tibble(long = c(-71.058884), lat = c(42.360081), names = c("MA\n82"))
  OrlandoBrew <- tibble(long = c(-81.379234), lat = c(28.538336), names = c("FL\n58"))
  NewYorkBrew <- tibble(long = c(-76.147423), lat = c(43.048122), names = c("NY\n74"))
  PortlandBrew <- tibble(long = c(-122.6750), lat = c(45.5051), names = c("OR\n125"))
  IndianaBrew <- tibble(long = c(-86.1581), lat = c(39.7684), names = c("IN\n139"))
  #View(map.df)
  
  # Plot the data on map
  p1 <- ggplot(map.df, aes(long, lat))+
    geom_polygon(aes(fill = BrewriesCnt, group = group), stat= "identity", color = "white")+
    #geom_text(aes(label=state), hjust=0.5, vjust=0.5, size=3)  
    #geom_text(aes(label = BrewriesMPD$count, )) +
    scale_fill_gradientn (colours= (heat.colors(3)),na.value="grey90")+
     ggtitle( "Breweries by State", subtitle = "Geo-positional Analysis") +
    coord_quickmap()
  p1 + geom_text(data = DenverBrew, aes(x = long, y = lat, label = names), hjust = 0.5, vjust=0.8,nudge_x = 0, color = "black", size=3) +
    geom_text(data = LosAngelesBrew, aes(x = long, y = lat, label = names), hjust = 0.5, vjust=-0.1,nudge_x = 0, color = "black", size=3) +
    geom_text(data = TraverseCityBrew, aes(x = long, y = lat, label = names), hjust = 0, vjust=0.8,nudge_x = 0, color = "black", size=3) +
    geom_text(data = AustinBrew, aes(x = long, y = lat, label = names), hjust = 0.8, vjust=-0.1, nudge_x = 0, color = "black", size=3) +
    geom_text(data = PhiladelphiaBrew, aes(x = long, y = lat, label = names), hjust = 1, vjust=-0.1,nudge_x = 0, color = "black", size=3) +
    geom_text(data = ChicagoBrew, aes(x = long, y = lat, label = names), hjust = 1, vjust=1, nudge_x = 0, color = "black", size=3) +   
    geom_text(data = BostonBrew, aes(x = long, y = lat, label = names), hjust = 0, nudge_x = 0, color = "black", size=3) +    
    geom_text(data = OrlandoBrew, aes(x = long, y = lat, label = names), hjust = 0.5, nudge_x = 0, color = "black", size=3) +
    geom_text(data = NewYorkBrew, aes(x = long, y = lat, label = names), hjust = 0.5, vjust=-0.2, nudge_x = 0, color = "black", size=3) +
    geom_text(data = PortlandBrew, aes(x = long, y = lat, label = names), hjust = 0, vjust=1, nudge_x = 0, color = "black", size=3) +
    geom_text(data = IndianaBrew, aes(x = long, y = lat, label = names), hjust = 0, nudge_x = 0, color = "black", size=3)

# Analyze External Datasource AlcoholDF
  AlcoholDF$region <- tolower(AlcoholDF$State)
# AlcoholDF$NumTot <- as.numeric(AlcoholDF$Total)
  AlcoholDF$NumTot <- as.numeric(AlcoholDF$GreaterThan12Est)
  AlcoholDF$chaTot<- as.character(AlcoholDF$NumTot)
  AlcoholMap.df <- left_join(states,AlcoholDF, by = "region")
  
  # Alcohol Abuse Top10 cities latitude and longitude
  LosAngeles <- tibble(long = c(-118.243683), lat = c(34.052235), names = c("CA\n7726"))
  Austin <- tibble(long = c(-97.743057), lat = c(30.267153), names = c("TX\n5302"))
  Orlando <- tibble(long = c(-81.379234), lat = c(28.538336), names = c("FL\n4370"))
  NewYork <- tibble(long = c(-76.147423), lat = c(43.048122), names = c("NY\n4201"))
  Chicago <- tibble(long = c(-87.629799), lat = c(41.878113), names = c("IL\n2949"))
  Philadelphia <- tibble(long = c(-75.165222), lat = c(39.952583), names = c("PA\n2909"))
  Columbus <- tibble(long = c(-82.998795), lat = c(39.961178), names = c("OH\n2425"))
  Boston <- tibble(long = c(-71.058884), lat = c(42.360081), names = c("MA\n1851"))
  Trenton <- tibble(long = c(-74.759720), lat = c(40.220581), names = c("NJ\n1836"))
  Atlanta <- tibble(long = c(-84.387985), lat = c(33.748997), names = c("GA\n1757"))
  Raleigh <- tibble(long = c(-78.638176), lat = c(35.779591), names = c("NC\n1728"))
  
  
  p2 <- ggplot(AlcoholMap.df, aes(long, lat)) +
    geom_polygon(aes(fill = GreaterThan12Est, group=group), stat= "identity", color = "white") +
    #geom_text(aes(label=state), hjust=0.5, vjust=0.5, size=3)  
    #geom_text(aes(label = BrewriesMPD$count, )) +
    scale_fill_gradientn (colours= (heat.colors(3)),na.value="grey90") +
    ggtitle( "Alcohol Abuse by State", subtitle = "Geo-positional Analysis") +
    coord_quickmap()
  p2 + geom_text(data = LosAngeles, aes(x = long, y = lat, label = names), hjust = 0.5, vjust=-0.1,nudge_x = 0, color = "black", size=3) +
    geom_text(data = Austin, aes(x = long, y = lat, label = names), hjust = 0.8, vjust=-0.1, nudge_x = 0, color = "black", size=3) +
    geom_text(data = Orlando, aes(x = long, y = lat, label = names), hjust = 0.5, nudge_x = 0, color = "black", size=3) +
    geom_text(data = NewYork, aes(x = long, y = lat, label = names), hjust = 0.5, vjust=-0.2, nudge_x = 0, color = "black", size=3) +
    geom_text(data = Chicago, aes(x = long, y = lat, label = names), hjust = 1, vjust=1, nudge_x = 0, color = "black", size=3) +
    geom_text(data = Philadelphia, aes(x = long, y = lat, label = names), hjust = 1, vjust=-0.1,nudge_x = 0, color = "black", size=3) +
    geom_text(data = Columbus, aes(x = long, y = lat, label = names), hjust = 0, nudge_x = 0, color = "black", size=3) +
    geom_text(data = Boston, aes(x = long, y = lat, label = names), hjust = 0, nudge_x = 0, color = "black", size=3) +
    geom_text(data = Trenton, aes(x = long, y = lat, label = names), hjust = 0, nudge_x = 0, color = "black", size=3) +
    geom_text(data = Atlanta, aes(x = long, y = lat, label = names), hjust = 0, nudge_x = 0, color = "black", size=3) +
    geom_text(data = Raleigh, aes(x = long, y = lat, label = names), hjust = 0, nudge_x = 0, color = "black", size=3)

# Merge Breweris and AlcoholAbuse dataset
  names(BrewriesMPD)[1] <- "State"
  BreweriesAlcoholAbuse <- left_join(BrewriesMPD, AlcoholDF, by = "State")

## Warning: Column `State` joining factors with different levels, coercing to
## character vector

  BreweriesAlcoholAbuse_new <- BreweriesAlcoholAbuse %>% select(State, BrewriesCnt, GreaterThan12Est)
  names(BreweriesAlcoholAbuse_new)[3] <- "AlcoholAbuseCnt"
  head(BreweriesAlcoholAbuse_new[order(-BreweriesAlcoholAbuse_new$AlcoholAbuseCnt),], n=11)

##            State BrewriesCnt AlcoholAbuseCnt
## 5     California         183            7726
## 43         Texas         130            5302
## 9        Florida          58            4370
## 32      New York          74            4201
## 13      Illinois          91            2949
## 38  Pennsylvania         100            2909
## 35          Ohio          49            2425
## 22      Michigan         162            2120
## 21 Massachusetts          82            1851
## 30    New Jersey           8            1836
## 10       Georgia          16            1757

  BreweriesAlcoholAbuse_new %>% ggplot(aes(x=BrewriesCnt, y=AlcoholAbuseCnt)) + geom_point() +
    ggtitle("Alcohol Abuse Count vs. Breweries Count in State") +
    geom_smooth(method = lm)

  cor(BreweriesAlcoholAbuse_new$BrewriesCnt, BreweriesAlcoholAbuse_new$AlcoholAbuseCnt)

## [1] 0.570742

Summary

Through our analysis:…